Directional interpolation for magnetic resonance angiography

ABSTRACT

The staircase artifact appearing in Magnetic Resonance Angiography boundaries between flowing blood and stationary tissue is eliminated by using a variable local interpolation direction which is determined based on comparing local intensity patterns. This technique of local directional interpolation is applied either to a three dimensional array of voxels prior to projection or to a two dimensional array of pixels after projection. In the case of interpolating lines of further pixels between consecutive lines of original pixels, a window of n 1  pixels in a present line is compared to successive groups of n 1  pixels in a search window of n 2  ≧n 1  pixels which groups are successively shifted in position in the search window. A match measure is formed by either a correlation or a root mean squared error (or a combination thereof) between the window and groups compared. The best match with respect to a threshold determines corresponding pixels. In the absence of a match better than the threshold a default vertical direction of interpolation is taken.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to angiography wherein a two dimensionalarray of pixels for display as an angiogram is determined from a threedimensional array of voxels computed from signal samples of radiationdue to flowing blood in a region of a body under examination. In itsparticular aspects, the present invention relates to interpolationmethods of providing an increased apparent resolution in the angiogram.

2. Background of the Invention

Display algorithms for rendering three-dimensional Magnetic ResonanceAngiography (MRA) data in two-dimensional form are known from H. Clineet al, "Volume Rendering and Connectivity Algorithms for MR Angiography"Magn. Res. Med. 18, pp. 384-394 (1991).

Angiography of the type mentioned is presently carried out by operatingmagnetic resonance imaging (MRI) apparatus in a mode employing amagnetic resonance angiography (MRA) technique in which a threedimensional data set of voxels exhibiting enhanced contrast to flowingblood is post-processed to display the vascular area of interest.Time-of-flight techniques for both 2D and 3D collections and 3Dphase-contrast techniques are known for enhancing the contrast offlowing blood relative to stationary tissue. In the 2D time-of-flightmethod, a collection of spin resonance signals for multiple parallelslices is obtained. The flow sensitive contrast is due to substantiallysaturating stationary spins in the slice from which spin resonancesignals are collected by relatively rapid repetition of close to 90°flip angle slice selective RF excitation pulses so that only unsaturatedspins in blood flowing into the slice have relatively stronglongitudinal magnetization just prior to the excitation pulses. Thisinduces high intensity spin resonance signals from the inflowing blood,which intensity increases with the amount of inflow velocity componentnormal to the slice. A three dimensional data set of voxel intensitiesis computed by two dimensional Fourier transformation for each slice ofsamples of the spin resonance signals received during a read gradientfor sequences repeated with different phase encoding gradient integrals.In the 3D time-of-flight method and in the 3D phase-contrast method athree dimensional data set of spin resonance samples are obtained fromwhich the three dimensional data set of voxel intensities is obtained bythree dimensional Fourier transformation.

Irrespective of the technique employed to obtain the three dimensionaldata set of voxel intensities, a rendering of this data set for viewingpurposes to a two dimensional data set of pixels is necessary. This isdone typically by forming a projection in a viewing direction. The mostwidely used projection method is maximum intensity projection (MIP).With this method, which is computationally fast, parallel rays areprojected through the three dimensional data set in a viewing direction,a different ray being associated with each pixel, and the maximumintensity of voxels along each each ray is taken as the intensity of theassociated pixel.

If the center to center spacing in one of three mutually orthogonaldirections, e.g. between adjoining voxels in different slices of amultislice collection, is larger than that between adjoining voxels inthe other directions, it is common to provide additional apparentresolution in the slice direction by interpolation. Interpolation isusually applied prior to projection increasing the number of planes ofvoxels in order to produce an expanded three dimensional data set withsubstantially cubic voxels. Alternatively, interpolation can be appliedto the two dimensional pixel array produced by projection to expand thenumber of lines quantizing the slice direction.

Both interpolation techniques produce an objectional staircase artifactalong the boundaries between flowing blood and stationary tissue whenblood vessels are oriented obliquely with respect to the principaldirections of the voxel spacings. This is because conventionalinterpolation techniques utilize a constant direction of interpolationwhich is aligned with one of the principle axes of the voxel array.

SUMMARY OF THE INVENTION

It is an object of the present invention to provide an angiographymethod using an interpolation technique which substantially minimizes oreliminates production of the staircase artifact. It is a further objectof the present invention to provide techniques for interpolation, eitherbefore or after projection, where local directions of interpolation aredetermined which tend to be along actual or projected longitudinal axesof oblique blood vessels. It is still another object of the inventionthat said local directions of interpolation be determined in acomputationally efficient manner not necessitating actual labelling ofpixels or voxels corresponding to flowing blood or vessel inner walls.

Briefly, the above objects are satisfied by providing an angiographymethod of the type where interpolation is done after projection of thethree dimensional array of voxels into a preliminary two dimensionalarray of pixels as a function of the positions and intensities of themeasured voxels. In such instance, additional lines of pixels areinterpolated between parallel lines of pixels in the preliminary twodimensional array in forming an expanded two dimensional pixel array fordisplay as an angiogram. In accordance with the invention, suchinterpolating is carried out by determining for each pixel in a firstline whether a pixel exists in an adjoining parallel said second linewhich corresponds in position to the pixel in the first line withrespect to intensity patterns in the lines, which are due to anyproximate cardiovascular structure, and by using a local direction ofinterpolation directed between the pixel in said first line and a pixelin the second line determined to correspond thereto, and in the absenceof a pixel in the second line determined to correspond thereto using alocal direction of interpolation along a perpendicular to said first andsecond lines.

Where the angiography method is such that interpolation is done prior toprojection, which may be preferable if a plurality of projections indifferent viewing directions are to be produced, then additional planesof voxels are interpolated between voxels of the three dimensional arrayto form an expanded three dimensional array to which the projectionalgorithm is applied. In this instance the determination ofcorrespondence relative to intensity patterns is between voxels inadjoining parallel planes, and the local interpolation direction, in thecase of correspondence, may have components in the three principaldirections of the voxel array.

The determination of correspondence is by forming a match measure whichis a function of the intensities of a first group of contiguous pixelsor voxels, in a first line or plane, including the pixel or voxel, in apredetermined relative position such as center, for which acorresponding pixel or voxel in the second line or plane is sought, andsecond groups of contiguous pixels or voxels, which second groups aresuccessively shifted in position in said second line or plane.Correspondence is determined with respect to pixel or voxel in the samepredetermined relative position in that second group, if any, whichyields the optimum match measure better than a predetermined thresholdvalue. The match measure is taken as a root mean squared difference ofintensities of the first and second groups, a correlation coefficient ofthe intensities and average intensities of the first and second groups,or a weighted combination of these alternatives.

Where no match measure is better that the predetermined threshold, thelocal direction of interpolation is taken, as in the prior art,perpendicular to the interpolated line of pixels or normal to theinterpolated plane of voxels.

Other objects, features and advantages of the present invention willbecome apparent upon study of the following detailed description whentaken in conjunction with the appended drawing wherein:

BRIEF DESCRIPTION OF THE DRAWING

FIG. 1 is a simplified schematic of magnetic resonance imaging apparatususable for implementing the present invention;

FIG. 2 is an illustration of the production of an angiogram byprojection;

FIG. 3 is a two dimensional representation of local directions ofinterpolation relative to blood vessels in an 2D array of pixelsobtained after projection, according to the principles of the presentinvention;

FIG. 4A is a three dimensional representation of oblique interpolationrelative to a 3D array of voxels prior to projection, in accordance withthe principles of the present invention;

FIG. 4B is a plan view of a plane in FIG. 4A;

FIG. 5 is a pictorial representation showing a first group of pixels orvoxels and a search window containing second groups of pixels or voxelswhich are compared to the first group in determining correspondence;

FIG. 6 is pictorial representation similar to FIG. 5, but with thecenter of the search window center biased in position; and

FIG. 7 is a flow chart of the directional interpolation method of thepresent invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

FIG. 1 of the drawing shows a simplified schematic of a magneticresonance imaging (MRI) apparatus 10 for the purpose of orientation.Therein, a main field coil 12 establishes a steady uniformlongitudinally directed magnetic field in an elongated measuring space14. The approximate range of field strengths in present MRI equipment isin the range of 0.1 to 4.0 Tesla. Main field coil 12 may besuperconducting, rather than resistive, in which case power supply 16 isused solely for initiating the necessary coil current. Also acting uponmeasurement space 14 are a system of three independently controllablegradient field coils 18 producing longitudinally directed fields withstrengths that varies linearly in respective three orthogonaldirections, in proportion to the respective coil currents driven bypower supply 20. An RF coil 22 receives RF pulses formed in RFtransmitter section 24 via Transmit/Receive switch 26 and transmitsthese pulses into measurement space 14. As is well known, these RFpulses are at the Larmor frequency determined by the strength of themain field and the nuclei of interest (usually water protons). RF coil22 receives the resultant spin resonance signals from the measurementsapce and via switch 26 feeds these signals to receiver section 28 wherethey are detected and then sampled and digitized by an A/D converter 30.A sequence controller controls gradient power supply 20,Transmit/Receive switch 26, RF transmitter section 24 and RF receiversection 28 so as to subject measurement space 14 to a sequence ofsub-experiments, each including at least a read or excitation RF pulsefor converting some or all longitudinal magnetization to transversemagnetization, and gradient pulses functioning as one or more phaseencoding gradients (the time integral of which is stepped from eachsub-experiment to the next) followed by a read gradient during whichspin resonance signals are sampled by the A/D converter 30.

Typically, for angiography using a time-of-flight procedure, thesequence is configured to induce spin resonance signals of relativelyhigh amplitude from voxels through which blood is flowing than fromvoxels containing only stationary tissue. A so-called 2D "multislice"collection is commonly utilized in which a slice-selective sequence isapplied sequentially to a stack of a plurality of (for example, four)equally spaced apart parallel slices, of known thickness. In aslice-selective sequence, at least the RF read pulse is appliedsimultaneously with a slice-selection gradient. The result of themultislice collection is a two-dimensional array (containing, forexample, 256×256 points in a spatial frequency space known as K-space)of signal samples for each slice, so that the multi-slice collection iseffectively a three dimensional array of signal samples which aresupplied by A/D converter 30 to digital memory 32. The signal sampledata for each slice is separately subjected to a two dimensional FourierTransform in block 34, which may constitute an array processor, toproduce a two-dimensional array of pixel intensities for each slice. Asa result, the collection of pixels for multiple parallel slices suppliedby block 34 to digital memory 36 can be considered a three-dimensionalarray of voxel intensities (of, for example, 4×256×256 points).

In the case of the use of either a 3D time-of-flight method or 3-Dphase-sensitive method, a three-dimensional array of sample points inK-space is read into memory 32 with usually one of the dimensions, whichis referred to herein as vertical, being quantized by less points thatthe other dimensions. This three-dimensional K-space array is thensubjected in block 34 to a three-dimensional Fourier Transform toproduce the three dimensional array of voxel intensities in memory 36.

This three-dimensional array of voxel intensities is operated upon by animage processor 38 to supply to digital memory 40 a two-dimensionalarray of pixel intensities which are displayed on a display device 42,such as a CRT. Image processor includes program modules for variousfunctions including background suppression 38a, interpolation 38b andprojection 38c, which are pertinent to angiography.

FIG. 2 illustrates application of projection module 38c to produce atwo-dimensional array of pixels for display as an angiogram 44 which isobtained by projecting parallel rays in a viewing direction 46 through athree-dimensional array 48 of voxels which are pixels of four stackedslices 50a through 50d, each ray producing an associated pixel of theangiogram. For purposes of illustration, a three dimensional intensitypattern corresponding to flowing blood 52 within a vessel is positionedintersecting three dimensional array 48. As a result of the projectionprocedure, the voxels making up the disc shaped volumes 54a, 54b ofvessel branches 52a, 52b within slice 50a, project respectively topixels making up the areas 54a, 54b of angiogram 44. The commonprojection technique for MR angiography is maximum intensity projection(MIP). In MIP, the intensity of each pixel in the angiogram is taken asthe maximum intensity of the voxels through which the associated raypasses. As, is apparent, in the absence of interpolation either beforeor after projection, the angiogram produced in a horizontal viewingdirection would be only four voxels high.

FIG. 3 illustrates a linear interpolation operation carried out byinterpolation module 38b when applied to interpolate, for example, apair of additional lines n1, n2 of pixels between each original pair ofadjoining parallel lines m, m+1 of pixels in a preliminary angiogramproduced by projection module 38a. As shown with respect to an obliquevessel branch 52a and a vertical vessel branch 52b, local interpolationdirections along lines 56 are utilized within and proximate to vessel orbranch 52. These local interpolation directions, rather than beingvertical, are oblique and tend to be parallel with the axis of branch52. These local directions are based on a determination ofcorrespondence of pixels in line m with pixels in line m+1, and theconnection of these corresponding pixels with lines 56. For example,pixels A and C in line m respectively correspond with pixels B and D inline m+1 based on their positions in the intensity patterns along thelines. By connecting the centers of pixels A and C respectively to thecenters of pixels B and D with straight lines 56, the points E, F ofintersection of lines 56 with line n1 are determined. Points E and F aretypically intermediate actual pixel centers. Then, hypotheticalintensities I_(E), I_(F) of points E, F respectively and therefrom, theintensity I_(G) of pixel G on line n1 between points E and F aredetermined from the following relationships, where I_(A) through I_(D)are the intensities of pixels A through D:

    I.sub.E =[I.sub.A ×(T-t)+I.sub.B ×t]/T

    I.sub.F =[I.sub.C ×(T-t)+I.sub.D ×t]/T

    I.sub.G =[I.sub.E ×L.sub.2 +I.sub.F ×L.sub.1 ]/(L.sub.1 +L.sub.2)

where as shown in FIG. 3:

t is the distance between line m and the line containing G;

T is the distance between lines m and m+1;

L1 is the distance between points E and G; and

L2 is the distance between points G and F.

Thus, with the oblique direction of interpolation illustrated, theintensity of any pixel G of an interpolated line is a function of theintensities of four pixels, A through D, of the original lines below andabove the interpolated line. It is this oblique direction ofinterpolation which tends to eliminate the staircase artifact present inthe edges of oblique vessels as they appear in MR angiograms produced byinterpolation.

It should be understood that the intensities I_(E) and I_(F) of theintermediate points E and F could be found from the corresponding pixelsin more than two lines, for example the corresponding points in a linem-1 and line m+2 (not shown) could be used at the same time and a thirdorder spline fit to the intensities at the corresponding points. In suchcase, the intensity I_(G) determined for an interpolated pixel G is afunction of the intensities of two pixels from each of the four originallines (m-1, m, m+1 and m+2) used to determine the intensities I_(E) andI_(F) of intermediate points E and F.

To interpolate pixels not within or proximate an oblique branch orvessel, a vertical interpolation procedure is used which is conventionalalthough contrary to the principles of the present invention, it washeretofore used for all pixels irrespective of their location relativeto anatomy or intensity patterns. By vertical interpolation, theintensity I_(L) of an interpolated pixel L is determined from theintensities I_(J) and I_(K) of the two original pixels J (line m) and K(line m+1) directly below and above pixel L, from the followingrelationship, where t and T have the same definitions as previouslygiven:

    I.sub.L =[I.sub.J ×(T-t)+I.sub.K ×t]/T

FIG. 4A illustrates interpolation of additional planes or slices n1 andn2 between original planes or slices m and m+1 of the three dimensionalarray of voxels prior to projection. It is an extension of thediscussion concerning FIG. 3 to three dimensions. An oblique directionof interpolation is shown where voxels D, E and F forming a triangle DEFin plane m+1 have been determined based on intensity patterns tocorrespond respectfully to voxels A, B and C forming a right triangleABC in plane m. The lines 56 between corresponding voxels of these twotriangles intersect plane n1 at points G, H and I forming a triangleGHI.

As better shown in FIG. 4B, points G, H and I are in generalintermediate actual voxel points to be interpolated. The hypotheticalintensities of points G, H and I and therefrom, the intensity of a voxelJ within triangle GHI are determined as an extension of therelationships discussed with respect to FIG. 3 as a function of thebetween plane distances t and T. One such extension is mathematicallyequivalent to averaging the intensities of the interpolated voxels byseparately applying the relationships discussed with respect to FIG. 3in each of the two in-plane orthogonal directions, using L₁, L₂ in onecalculation and L₃, and L₄ in the other. Alternatively, the intensity ofpixel J may be determined as a function of the distances M₁, M₂ and M₃from voxels J to points G, H and I.

As will become clearer as the discussion proceeds, a key aspect of thepresent invention is the determination of corresponding pixels indifferent lines or corresponding voxels in different planes based onlocal intensity patterns. This correspondence is determined by a searchprocess illustrated in FIG. 5, which is conveniently discussed in termsof determination of corresponding pixels in adjoining lines in a twodimensional array applicable to interpolation after projection. Itshould be understood that this discussion is equally applicable to thedetermination of corresponding voxels in adjoining planes in a threedimensional array, applicable to interpolation prior to projection. Theextension of a one dimensional search process to two dimensions is selfevident.

FIG. 5 illustrates the matching process between two lines, where thelines are denoted by m and m+1. The goal is to find for each point online m the corresponding point on line m+1. The figure illustratesmatching for the point at x coordinate x₀. This is performed by placinga reference window of size n₁ that is centered at the point x₀, on linem. A larger search window of size n₂ is centered at the point x₀, online m+1. These windows are represented by P_(m) (X) and P_(m+1) (X),respectively. The first number in parenthesis denotes the centercoordinate of the window. A number of measures can be used to estimatethe degree of the match. For a pair of windows the root mean squareerror (rmse) or the correlation coefficient can be calculated for everyshift value, s, where the smaller window has a complete overlap with thelarger window. The shift, s, is relative to the center of the referencewindow. The rmse is then calculated by: ##EQU1## where the range for thewindow shifts is: ##EQU2## and the normalized rmse is given by: ##EQU3##where e(s) ε [0,1]. The best match is the value of s for which e isminimum.

Similar to the rmse, the correlation coefficient can be calculated atevery point of complete overlap between the two windows. Let P_(m) andP_(m+1) represent the average values of pixels in the two windows. Let:

    W.sub.m (X)=P.sub.m (X)-P.sub.m (X)

    and:

    W.sub.m+1 (x)=P.sub.m+1 (x)-P.sub.m+1 (x)

The correlation coefficient is then given by: ##EQU4##

The correlation coefficient ρ is bounded in the range [-1,1]. It is +1for positively correlated and -1 for negatively correlated windows. Thecorrelation coefficient is insensitive to window intensity averages.Thus, a transition from low to mid-gray can be matched to one frommid-gray to bright. The normalized rmse is sensitive to average windowintensity levels, but is insensitive to changes in contrast. Acombination of correlation coefficient and normalized rmse may provide amore robust method for template matching. Combined measures can be(ρ-e), or weighted differences of the two measures, etc.

The matching and interpolation method, as described, can cause adirectional texture appearance in the background noise regions. Acontrol parameter in the algorithm was designed to avoid this problem.The background noise in MRA images is random and can be approximated bya Gaussian. The strength of the match in a background region is oftenmuch weaker than in regions containing flow. Therefore, one can avoidartifacts by looking at the degree of the match and comparing it with apreset threshold. If the match is better than the threshold it isaccepted. Otherwise, the calculated match is ignored and the defaultvertical match is assumed. For example, the normalized rmse can lie inthe range 0 to 1, where the strongest match is 0, and the weakest is 1.Suppose that normalized rmse is the match criterion, and that athreshold of 0.5 is chosen. Then all matches with rmse values greaterthan 0.5 are assumed to be in background regions, and a zero shiftvertical match is imposed. If the underlying assumption that matchstrengths in flow regions are better than the threshold is violated,directional interpolation is not carried out for those flow regions aswell.

A variation on the method is to use the match information from previouslines to bias the direction of the search. As shown in FIG. 6 is anapproach m+1 to points on line m+2. In the previous line matchingiteration, the point x₀ on line m was matched to the point x₀₋₂ on linem+1. This match represents a shift Δ=-2. The center of the search windowon line m+2 is shifted by the value of the bias, Δ, which is -2 in thisfigure. The search window is, therefore, centered on the point x₀₋₄. Inthe general case of such an approach the range for possible shifts, s,is given by: ##EQU5##

The bias effectively increases the search range in one direction, whilereducing it in the other. Given a fixed search window size, the biasincreases the change of finding the right match if the vessel shiftslarger distances in the same direction as the previous line. However, itreduces the chance of finding the right match if the vessel suddenlycurves in the opposite direction. An alternative (not shown in FIG. 7)is to adaptively increase or decrease the size of the search windowbased on prior matches. For example, if the previous match is close tothe edge of the search window, the search window size is increased.Another variation is to use the matching information from the previouspoint on the same line to bias the match.

FIG. 7 is a flow chart, which is self explanatory, of the method ofdetermining corresponding pixels which determines one or moreinterpolated lines between present and next lines of a projection imageonce all corresponding pixels have been found in these present and nextlines. As shown, first in block 60 window sizes and threshold for thematch measure are selected. The window sizes thus far found preferableare n₁ =11 and n₂ =17. However these parameters should be optimized forparticular anatomy and protocols. An appropriate threshold for the matchmeasure depends on whether background suppression is performed prior tointerpolation and should also be optimized for the particular anatomyand protocol.

In blocks 62 and 64, a first and sequential next lines of pixels areread. In block 66, there is a move to the first allowed point on thepresent line for reference matching. Blocks 68 through 80 are within aninner loop which sequentially determines corresponding points in thenext line to the present line. In block 70 through 74, the search may becarried out either biased as illustrated in FIG. 6 or not biased as inFIG. 5, depending on the selection made in box 70. The optimum matchmeasure selected in box 76 is preferably the closest of multiplethereof. In block 78, when a match is rejected, the default of avertical match is taken. After the inner loop comprising blocks 68through 80, an interpolated line is determined in block 84. As shouldnow be apparent, an outer loop comprising blocks 64 through 88 inresponse to each read in of a next line, determines the one or moreinterpolated lines between the present line and this next line.

As previously noted, it should be apparent that the flow chart of FIG. 7is easily extended to three dimensions by any desired pattern ofsearching over the shifts in each of two dimensions for the optimummatch measure.

It should also be apparent that although the present invention has beendescribed in detail with respect to MR angiography, many principlesthereof are applicable to other medical imaging modalities such as CT.Further, numerous modifications in the details given are possible withinthe intended spirit and scope of the invention.

What is claimed is:
 1. An angiography method comprising:inducingradiation to radiate from voxels of a region of a body under examinationhaving intensities influenced by any flowing blood within said voxels,said flowing blood being in blood vessels, said voxels having center tocenter spacing in respective three orthogonal directions; receiving andsampling the radiation radiating from the voxels of said region as acollection of signal samples; converting said collection of signalsamples into a three dimensional array of computed voxel intensities,each dimension of said array corresponding to a different one of saidthree orthogonal directions; forming a preliminary two dimensional arrayof pixels as a function of the positions and intensities of the measuredvoxels; and interpolating between first and second parallel lines ofpixels in the preliminary two dimensional array, an additional line ofinterpolated pixels in forming an expanded two dimensional pixel arrayfor display as an angiogram; wherein said interpolating is carried outby determining for each pixel in said first line whether a pixel existsin said second line which corresponds in position to the pixel in thefirst line with respect to local intensity patterns, and by using alocal direction of interpolation directed between the pixel in saidfirst line and a pixel in the second line determined to correspondthereto, and in the absence of a pixel in the second line determined tocorrespond thereto using a local direction of interpolation along aperpendicular to said first and second lines; whereby local directionsof interpolation tend to be parallel to axes of proximate blood vessels.2. The method as claimed in claim 1, wherein said determining comprisesforming a match measure between a first group of consecutive pixels insaid first line, which first group includes in a predetermined relativeposition in said first group, the pixel for which the existence of acorresponding pixel in said second line is being determined, and each ofa plurality of second groups of consecutive pixels, which second groupsare successively shifted in position in said second line, andidentifying as a corresponding pixel in said second line, the pixel inthe same predetermined relative position in the second group having abest match measure of the plurality of second groups, if said best matchmeasure is better than a predetermined threshold value.
 3. The method asclaimed in claim 2, wherein said forming of said match measure withrespect to a second group comprises forming a sum over the respectivesame relative positions in the first and second groups, of the square ofthe differences between intensities of pixels having the same relativepositions in the respective first and second groups.
 4. The method asclaimed in claim 3, wherein said forming of said match measure withrespect to a second group comprises forming a sum over the respectivesame relative positions in the first and second groups, of products ofintensities of pixels at the same relative positions in the respectivefirst and second groups biased by mean intensities applicable to saidrespective first and second groups.
 5. The method as claimed in claim 2,wherein said forming of said match measure with respect to a secondgroup comprises forming a sum over the respective same relativepositions in the first and second groups, of products of intensities ofpixels at the same relative positions in the respective first and secondgroups biased by mean intensities applicable to said respective firstand second groups.
 6. The method as claimed in claim 2, wherein saidforming a preliminary two dimensional display array of pixels is bydetermining the maximum intensity of voxels along parallel rays directedin a viewing direction;
 7. The method of claim 2, wherein said inducedand received radiation are MR spin resonance RF signals.
 8. The methodas claimed in claim 1, wherein said forming a preliminary twodimensional display array of pixels is by determining the maximumintensity of voxels along parallel rays directed in a viewing direction;9. The method of claim 1, wherein said induced and received radiationare MR spin resonance RF signals.
 10. An angiography methodcomprising:inducing radiation to radiate from voxels of a region of abody under examination with intensities influenced by any presence offlowing blood within said voxels, said flowing blood being in bloodvessels, said voxels having center to center spacing in respective threeorthogonal directions;; receiving and sampling the radiation radiatingfrom said voxels as a collection of signal samples; converting saidcollection of signal samples into a three dimensional array of computedvoxel radiation intensities, each dimension of said array correspondingto a different one of said three orthogonal directions; interpolating anadditional plane of voxels between first and second adjoining parallelplanes of said computed voxels in forming an expanded voxel array; andforming a two dimensional display array of pixels, for display as anangiogram, as a function of the positions and intensities of the voxelsin the expanded voxel array; wherein said interpolating is carried outby determining for each pixel in said first plane whether a pixel existsin said second plane which corresponds in position to the pixel in thefirst line with respect to local intensity patterns, and by using alocal direction of interpolation directed between the pixel in saidfirst plane and a pixel in the second plane determined to correspondthereto, and in the absence of a pixel in the second plane determined tocorrespond thereto using a local direction of interpolation along anormal to said first and second planes; whereby local directions ofinterpolation tend to align parallel to axes of proximate blood vessels.11. The method as claimed in claim 10, wherein said determiningcomprises forming a match measure between a first group of consecutivepixels in said first plane, which first group includes in apredetermined relative position in said first plane, the pixel for whichthe existence of a corresponding pixel in said second plane is beingdetermined, and each of a plurality of second groups of consecutivepixels, which second groups are successively shifted in position in saidsecond plane, and identifying as a corresponding pixel in said secondplane, the pixel in the same predetermined relative position in thesecond group having a best match measure of the plurality of secondgroups, if said best match measure is better than a predeterminedthreshold value.
 12. The method as claimed in claim 11, wherein saidforming of said match measure with respect to a second group comprisesforming a sum over the respective same relative positions in the firstand second groups, of the square of the differences between intensitiesof pixels having the same relative positions in the respective first andsecond groups.
 13. The method as claimed in claim 12, wherein saidforming of said match measure with respect to a second group comprisesforming a sum over the respective same relative positions in the firstand second groups, of products of intensities of pixels at the samerelative positions in the respective first and second groups biased bymean intensities applicable to said respective first and second groups.14. The method as claimed in claim 11, wherein said forming of saidmatch measure with respect to a second group comprises forming a sumover the respective same relative positions in the first and secondgroups, of products of intensities of pixels at the same relativepositions in the respective first and second groups biased by meanintensities applicable to said respective first and second groups. 15.The method as claimed in claim 11, wherein said forming a twodimensional array of pixels is by determining the maximum intensity ofvoxels along parallel rays directed in a viewing direction.
 16. Themethod of claim 11, wherein said induced and received radiation are MRspin resonance RF signals.
 17. The method as claimed in claim 10,wherein said forming a two dimensional array of pixels is by determiningthe maximum intensity of voxels along parallel rays directed in aviewing direction.
 18. The method of claim 10, wherein said induced andreceived radiation are MR spin resonance RF signals.
 19. Magneticresonance angiography apparatus comprising:means for inducing magneticresonance signals to radiate from voxels of a region of a body underexamination having intensities influenced by any flowing blood withinsaid voxels, said voxels having center to center spacing in respectivethree orthogonal directions; means for receiving and sampling themagnetic resonance signals radiating from the voxels of said region as acollection of signal samples; means for converting said collection ofsignal samples into a three dimensional array of computed voxelintensities, each dimension of said array corresponding to a differentone of said three orthogonal directions; means for forming a preliminarytwo dimensional array of pixels as a function of the positions andintensities of the measured voxels; and means for interpolating betweenfirst and second parallel lines of pixels in the preliminary twodimensional array, an additional line of interpolated pixels in formingan expanded two dimensional pixel array for display as an angiogram;wherein said interpolating means comprises means for determining foreach pixel in said first line whether a pixel exists in said second linewhich corresponds in position to the pixel in the first line withrespect to local intensity patterns, and for using a local direction ofinterpolation directed between the pixel in said first line and a pixelin the second line determined to correspond thereto, and in the absenceof a pixel in the second line determined to correspond thereto for usinga local direction of interpolation along a perpendicular to said firstand second lines; whereby local directions of interpolation tend to beparallel to axes of proximate blood vessels.
 20. MR angiographyapparatus comprising:means for inducing MR spin resonance RF signals toradiate from voxels of a region of a body under examination withintensities influenced by any presence of flowing blood within saidvoxels, said voxels having center to center spacing in respective threeorthogonal directions; means for receiving and sampling the MR spinresonance RF signals radiating from said voxels as a collection ofsignal samples; means for converting said collection of signal samplesinto a three dimensional array of computed voxel radiation intensities,each dimension of said array corresponding to a different one of saidthree orthogonal directions; means for interpolating an additional planeof voxels between first and second adjoining parallel planes of saidcomputed voxels in forming an expanded voxel array; and means forforming a two dimensional display array of pixels, for display as anangiogram, as a function of the positions and intensities of the voxelsin the expanded voxel array; wherein said interpolating means comprisesmeans for determining for each pixel in said first plane whether a pixelexists in said second plane which corresponds in position to the pixelin the first line with respect to local intensity patterns, and forusing a local direction of interpolation directed between the pixel insaid first plane and a pixel in the second plane determined tocorrespond thereto, and in the absence of a pixel in the second planedetermined to correspond thereto for using a local direction ofinterpolation along a normal to said first and second planes; wherebylocal directions of interpolation tend to be parallel to axes ofproximate blood vessels.